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Abstract 

We present a new method for reconstructing two-dimensional mass maps of galaxy clusters 
from the image distortion of background galaxies. In contrast to most previous approaches, 
' which directly convert locally averaged image ellipticities to mass maps (direct methods), our 

entropy-regularized maximum-likelihood method is an inverse approach. Albeit somewhat 
more expensive computationally, our method allows high spatial resolution in those parts of 
the cluster where the lensing signal is strong enough. Furthermore, it allows to straightfor- 
wardly incorporate additional constraints, such as magnification information or strong- lensing 
features. Using synthetic data, we compare our new approach to direct methods and find in- 
deed a substantial improvement especially in the reconstruction of mass peaks. The main 
differences to previously published inverse methods are discussed. 



o 

oo 

On 



q 1 Introduction 

The reconstruction of projected cluster mass maps from the observable image distortion of faint 
background galaxies due to the tidal gravitational field is a new and powerful technique. Pioneered 
by Kaiser & Squires (1993), this method has since been modified and generalized to account for 
(a) strong tidal fields in cluster centers (Schneider & Seitz 1995; Seitz & Schneider 1995; Kaiser 
1995); (b) finite and - in some cases, e.g. WFPC2 images - very small data fields (Schneider 1995; 
Kaiser et al. 1995; Bartelmann 1995; Seitz & Schneider 1996, 1998; Lombardi & Bertin 1998); and 
(c) the broad redshift distribution of background galaxies (Seitz & Schneider 1997). All of these 
are direct methods in the sense that a local estimate of the tidal field is derived from observed 
galaxy ellipticities, which is then inserted into an inversion equation to obtain an estimate of the 
surface mass density of the cluster. 

Whereas these direct methods are computationally fast, can be treated as black-box routines, 
need only the observed ellipticities and a smoothing length 8q as input data, and yield fair estimates 
of the surface mass density, their application has several drawbacks: 

• The data must be smoothed, and the smoothing scale is typically a free input parameter 
specified prior to the mass reconstruction. There arc no objective criteria on how to set 
the smoothing scale, although some ad-hoc prescriptions for adapting it to the strength of 
the lensing signal have been given (Seitz et al. 1996). In general, smoothing leads to an 
underestimate of the surface mass density in cluster centers or sub-condensations. 

• The quality of the reconstruction is hard to quantify. 

• Constraints on the mass distribution from additional observables (such as multiple images 
or giant arcs) cannot simultaneously be included. In particular, magnification information 
contained in the number density of background sources (Broadhurst et al. 1995; Fort et al. 
1997) or in the image sizes at fixed surface brightness (Bartelmann & Narayan 1995), cannot 
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be incorporated locally but only globally to break the mass-sheet degeneracy (Gorenstein et 
al. 1988; Schneider & Seitz 1995). 

To overcome these drawbacks, a different class of methods should be used. Bartclmann et al. 
(1996, hereafter BNSS) developed a maximum-likelihood (ML) technique in which the values of 
the deflection potential at grid points are considered as free parameters. After averaging image 
ellipticities and sizes over grid cells, local estimates of shear and magnification are obtained. The 
deflection potential at the grid points is then determined such as to optimally reproduce the 
observed shear and magnification estimates. Magnification information can be included this way. 
The smoothing scale in this method is given by the size of the grid cells, and can be chosen such 
that the overall % 2 of the fit is of order unity per degree of freedom. 

Squires & Kaiser (1996; hereafter SK) suggested several inverse methods. Their maximum 
probability method parameterizes the mass distribution of the cluster by a set of Fourier modes. 
If the number of degrees of freedom (here the number of Fourier modes) is large, the mass model 
tends to over-fit the data. This has to be avoided by regularizing the model, for which purpose 
SK impose a condition on the power spectrum of the Fourier modes. SK's maximum-likelihood 
method specifies the surface mass density on a grid and uses the Tikhonov-Miller regularization 
(Press et al. 1992, Sect. 18.5). The smoothness of the mass reconstructions can be changed by 
varying the regularization parameter, which is chosen such as to give an overall x 2 ~ 1 per degree 
of freedom. 

Bridle et al. (1998) have recently proposed an entropy-regularized ML method in which the 
cluster mass map is parameterized by the surface mass density at grid points. This method allows 
to restrict the possible mass maps to such with non-negative surface mass density. 

This paper describes another variant of the ML method (Seitz 1997, Ph.D. thesis). The major 
differences to the previously mentioned inverse methods are the following: 

• The observational data (e.g. the image ellipticities) are not smoothed, but each individual 
ellipticity of a background galaxy is used in the likelihood function. Whereas this modifica- 
tion complicates the implementation of the method, it allows larger spatial resolution for a 
given number of grid points, which is useful since the latter determines the computing time. 

• The number of grid points can be much larger than in BNSS, and the likelihood function 
is regularized. This produces mass reconstructions of variable smoothness: Mass maps are 
smooth where the data do not demand structure, but show sharp peaks where required by 
the data. The resulting spatially varying smoothing scale is a very desirable feature. Fourier 
methods, such as SK's maximum probability method, have a spatially constant smoothing 
scale which is determined by the highest-order Fourier components. They always need to 
compromise between providing sufficient resolution near mass peaks and avoiding over-fitting 
of the data in the outer parts of a cluster. 

• Following BNSS, we use the deflection potential to describe a cluster. This is an essential 
difference to Bridle et al. (1998) who used the surface mass density at grid points. As we 
shall discuss below, working with the deflection potential has substantial fundamental and 
practical advantages. 

We describe our method in Sect. 2, with details given in the Appendix. We then apply the 
method to synthetic data sets in Sects. 3 & 4 to demonstrate its accuracy. In particular, we compare 
the performance of our ML method to that of direct methods. The results are then discussed in 
Sect. 4, and conclusions are given in Sect. 5, where we also discuss further generalizations of the 
method for, e.g., including constraints from strong lensing features. 

2 The entropy-regularized ML mass reconstruction 
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2.1 Basic lensing relations 



For simplicity, we assume throughout the paper that all background galaxies are located at the 
same redshift. A generalization of our technique to a redshift distribution is given by Geiger & 
Schneider (1998). The dimensionless surface mass density n(8) is related to the deflection potential 
i/j(0) through the Poisson equation, 



(1) 



where indices i preceded by a comma denote partial derivatives with respect to 9i. The tidal 
gravitational field of the lens is described by the shear 7 = 71 + i 72 with the two components 



(2) 



Thus, the surface mass density and the shear, which determine the local properties of the lens 
mapping, can locally be obtained from the deflection potential ip. In contrast, the relation between 
shear and surface mass density is highly non-local, 



7(0) = - ( dV V{6 - 6') k{6') , 



(3) 



with V(6) = -{6\ - 0f + 2i6 1 6 2 )\9\- 4: . In particular, k needs to be given on the entire two- 
dimensional plane. Prescribing k on a finite field does therefore not completely specify the shear 
inside the field, because the latter is also affected by the outside mass distribution. We return to 
this point further below. The local magnification is 



(1(6)= { 1- K (6)] 2 - 7(0) 12 



(4) 



The local lens equation relates the ellipticities of a source and its image. We use the complex 
ellipticity parameter \ ( see Blandford et al. 1991) to describe image shapes. It is generally defined 
in terms of the tensor Qij of second brightness moments of an image by x = (Q11 ~ Q22 + 

2iQi2)/(Qii + Q 22 ). 



2.2 The ML method 

We refer the reader to Press et al. (1992, Chap. 18) for the basic ingredients of the ML method; 
see also Bridle et al. (1998). We do not repeat the basics here, but describe the application of 
the method to the cluster mass reconstruction. We start by considering image ellipticities only; 
magnification effects will be discussed later. 

Let Xk, 1 < k < N g , denote the complex ellipticities of N g galaxy images in the data field U, 
which we assume to be a rectangle of side lengths L x and L y . We cover the data field with an 
equidistant grid of N x x N y points &ij , with 6u in the lower left corner of the data field. The cluster 
is described by the deflection potential at the grid points, ipij, < i < N x + 1, < j < N y + 1. As 
discussed in BNSS, the grid for t/j is larger than the data field by one column or row of grid points 
in all four directions to allow simple finite differencing of ip on the whole field U. Having found k 
and 7 on all grid points from -0 according to ((!]) and ^), k and 7 are bilinear ly interpolated to all 
galaxy positions 8^ . 

If the isotropic probability distribution p s {x s ) °f the intrinsic source ellipticities is given, the 
probability distribution p(x'i {V'}) of t ne image ellipticities can be predicted. The likelihood func- 
tion is then defined as 

C(m):=l[p(xk;W). (5) 

fc=i 
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C, or ln£, can be maximized with respect to the set {ip} of values of the deflection potential at 
the grid points. Since the values of k and 7 are obtained from second derivatives of ip, a constant 
and a term linear in 9 added to ip leave C unchanged. In addition, the mass-sheet degeneracy 
renders C invariant under the transformation 

^9)^\m + ^^-\9\\ (6) 

where A ^ is an arbitrary parameter (Schneider & Seitz 1995). Therefore, in maximizing C with 
respect to {ip}, the potential tp can be held fixed at four grid points. Noting that the corners 
of the ip grid are not used for the calculation of n and 7 on U (see Appendix), we see that the 
maximization of C has dimension iV dim = (N x + 2)(N y + 2) — 8. 



2.3 Regularization 

Provided -/Vdim is not much smaller than the number of galaxies (which we assume in the following), 
the maximization results in a cluster model which tries to follow closely the noise pattern of the 
data. Disregarding observational effects, the noise is due to the intrinsic ellipticity distribution 
of the sources. The reconstructed mass distribution will therefore have pronounced small-scale 
structure, fitting the observed image cllipticities as closely as possible, and having a % 2 per degree 
of freedom much smaller than unity. In order to prevent such over-fitting of the data, we need to 
augment £ by a regularization term. Instead of maximizing ln£, we minimize 

E(M):=-^ln£(M)+»7ft(M), (7) 

where TZ is a function of the potential that disfavors strong small-scale fluctuations. The parameter 
rj determines how much weight should be attached to smoothness. One can vary r\ such that the 
resulting reconstruction has approximately the expected deviation from the data, viz. x 2 ~l per 
degree of freedom. Larger values of r\ yield mass distributions which are too smooth to fit the 
data, lower values of r\ cause over-fitting. 

We experimented with quite a large number of regularization terms. For example, we chose 
TZ as the sum of |Vk| 2 over all grid points. Mass reconstructions from synthetic data (see Sect. 3 
below) then showed a strong tendency to decrease too slowly towards the outer parts of the cluster, 
for that regularization preferred k to be as flat as possible. Rcgularizations including higher-order 
derivatives of k (see Press et al. 1992, Sect. 18.5) led to similar artifacts. Thus, such local linear 
rcgularizations were dismissed as unsatisfactory. 

Motivated by the success of the maximum-entropy (ME) image deconvolution (e.g. Narayan & 
Nityananda 1986; Lucy 1994), we consider instead ME regularizations of the form 

N x ,N y 

where 



Kij — 



N X ,N, 

k,i=i 



-1 -1 

Kij (9) 



is the normalized surface mass density at the grid points, and bij is a similarly normalized prior 
distribution (see Press et al. 1992, Sect. 18.7, for a detailed discussion of the ME method). We 
experimented with different choices for the prior. On the whole, a uniform prior, bij = const., 
performed satisfactorily, but tended to smooth mass peaks more than desired. Following Lucy 
(1994), we therefore use a prior which is determined by the data itself. Deferring details to 
Sect. 3, we note here that one can use the mass distribution obtained from a direct (finite field) 
reconstruction method as an initial prior, then iteratively minimize E, and after several iterations 
use a smoothed version of the current mass distribution as a new prior. Lucy (1994) showed that 
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such a moving prior yields more accurate reconstructions than a constant prior. The regularization 
parameter rj can be iteratively adjusted to provide the expected goodness-of-fit. Note that the 
ME regularization ensures that the reconstructed surface mass distribution is positive definite. 



2.4 The ellipticity distribution 

Given the intrinsic ellipticity distribution and the local values of k and 7, the probability distribu- 
tion p(x) for the image ellipticities can be calculated. However, the resulting analytic expressions 
are quite cumbersome and unsuitable for the high-dimensional minimization problem considered 
here. With the intrinsic distribution not being accurately known anyway, a precise expression of 
p(x) is not needed. We therefore approximate the image ellipticity distribution by a Gaussian, 
with mean (x) and dispersion a x . Both these values depend on the (local) distortion, 

5 = 29 - , (10) 

1 + M 2 

where g — 7/(1 — k) is the reduced shear. Mean and dispersion can be approximated by (Schneider 
& Seitz 1995) 



(x) = CS 



M9 / o\ Ml 

(i-i*r 



a x = <tq (l - |<5| 2 )" 2 , (11) 

where 

. 3M 4 6M 2 +Mf-9M 4 

Ml = 1 -4M^' ^ 2 = 8M 2 ' (12) 

and M n is the n-th moment of the intrinsic ellipticity distribution. Using (JSJ) and (Q), the function 
E to be minimized can then be written 

E(m)=G x m)+j2H4(k)]+vnm) , (is) 
is=i 

where we have introduced 

n — 1 \xk - (x) (fc)l 2 ri4 x 

x "^h ^ ' 

(x) (k) and <7 x (k) are the values of (x) and a x at the position of the fc-th galaxy, which obviously 
depend on the deflection potential {-0} through the distortion <5.[] The term G x in E has the 
form of a x 2 -function, which implies that an acceptable mass model should have G x ~ 1. This 
condition constrains the regularization parameter r\. 

We outline in the Appendix an efficient method for calculating E and its derivatives with 
respect to the values of ip 011 the grid points. We note that the terms in the sums of ( |l3| ) and 
( |l4| ) corresponding to a galaxy at 9k depend only on the values of tp at neighboring grid points. 
In that sense, our method is local. Had we parameterized the cluster by k at grid points, each 
term in the sums depended on k at all grid points, as can be seen from (||). It becomes obvious 
that the description in terms of ip requires much less computer time. In addition, the use of k 
rather than ip is strongly disfavored by the fact that the shear on U is incompletely specified by 
k on the same field. Bridle et al. (1998) attack this problem by performing the reconstruction on 
a region much larger than IA. They find that the shear information within 11 yields information 
about k outside hi. Although this is true, the information on the mass outside the data field 

1 If we worked in terms of the ellipticity parameter e (see BBNS) rather than the ellipticity \> then (e) = g (Seitz 
& Schneider 1997) for a non-critical cluster. In the general case including critical clusters, the variable \ ls more 
convenient. 
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is very limited: For a circular data field, the shear inside U caused by mass outside U can be 
fully described by conveniently defined multipole moments of the mass distribution outside hi 
(Schneider & Bartelmann 1997), and there are infinitely many mass distributions for which all 
of those multipole moments agree. For instance, a point mass located just outside the data field 
produces the same shear pattern in the field as a spherically symmetric mass distribution with the 
same total mass. Finally, extending the region on which the reconstruction is performed increases 
the dimensionality of the minimization problem. 

2.5 Including magnification information 

The mass-sheet degeneracy @ can be lifted if the lens magnification can be estimated. Three 
different methods for measuring magnification were suggested in the literature. Broadhurst et 
al. (1995) proposed to use the magnification bias, which changes the local number density of 
background galaxies due to the magnification, provided the slope of the number counts d log N/ dm 
is sufficiently different from 0.4. Noting that lensing magnifies objects but leaves their surface 
brightness unchanged, Bartelmann & Narayan (1995) suggested that the sizes of background 
galaxies at fixed surface brightness could be a convenient measure for the magnification after 
calibrating with field galaxies. Both methods can be used locally or globally. In the first case, 
the local magnification information is used for the mass reconstruction, whereas in the latter case, 
the transformation parameter A in (P) is adjusted until the magnification optimally matches the 
observational estimate. Kolatt & Bartelmann (1998) suggested to calibrate A globally by using 
type-la supernovae as cosmological standard candles. If magnification information is taken into 
account, the potential ip can be kept fixed at three points only, yielding iVdi m = {N x + 2)(N y +2) — 7. 

As an example, consider the method suggested by Bartelmann & Narayan (1995). Let r s and 
r be the ratios of the linear sizes of a galaxy and its image, respectively, relative to the mean size 

1/2 

of galaxies with the same surface brightness. They are related by r — \fi\ r s . The expectation 

1/2 

values of r s and r are 1 and \fjt\ , respectively. Hence, if p s (r s ) is the probability distribution for 
r s , the local probability distribution for r is 



Including the size distribution into the likelihood maximization leads to an additional term in 




(15) 



2?(M) 



G X (M) + ~ 2 Gr (M) 



+ 




(16) 



k=l 



where 




k=l 



(17) 



We assume p s to be a log-normal distribution (cf. Bartelmann & Narayan 1995), 




(18) 



with (r s ) = 1, and ([Tt]) becomes 




(19) 



G 



with (lnr) (k) = (In |/i(fc)| — of )/2. Hence, p is a Gaussian in lnr, with dispersion ay. G r has the 
form of a x 2 function in lnr, which motivates us to include the factor 1/2 in the definition of G r . 
A satisfactory mass model should have G r « 1. Note that the galaxies whose size can reliably be 
measured need not be those used for measuring the shear. It is just for notational simplicity that 
we assume that the two galaxy populations agree. 

3 Practical implementation and simulation parameters 

3.1 Practical implementation 

For minimizing ( |l3| ) or ( |l6| ) in the process of the ML mass reconstruction, the various quantities 
(e.g. (x) (k), <J x (k), kij) need to be calculated for each set of values {ip} of the deflection potential. 
We outline in the Appendix how this can be done efficiently. In order to quickly approach the 
minimum of E, we use derivative information in the minimization procedure. The derivative of 
E with respect to tpij is also given in the Appendix. We employ the conjugate gradient method 
as encoded in the routine frprmn by Press et al. (1992), with line minimization using derivative 
information. 

We need a good initial potential to start the minimization. We tested two different approaches. 
The first starts with a relatively small grid (say, N x = N y = 11 for a quadratic data field U), and 
a potential which corresponds to a constant surface mass density of k ~ 0.1, say. rj can be set to 
zero initially because the large grid cells provide sufficient smoothing to avoid over-fitting the data, 
or otherwise set 77 to a small value with a constant prior. After several (20, say) iterations, the 
current mass map is smoothed and used as the new prior, and the minimization is continued. The 
regularization parameter is slightly increased or decreased, depending on whether G x is smaller 
or larger than unity. Once a stable minimum with G x ~ 1 is obtained, the solution can be 
interpolated to a finer grid [N x = N y = 21, say) and the minimization can continue, adapting 
the prior and the regularization parameter as described before. This procedure can be repeated if 
desired. 

The second approach starts on a fine grid right away, with initial conditions obtained from a 
direct finite-field reconstruction method like that described by Seitz & Schneider (1998). From 
this mass distribution, an approximate deflection potential can be found by integration, although 
any contribution from the mass outside the data field will be missing in the resulting tp. This 
initial prior is a smoothed version of the mass map from the direct reconstruction, and is adapted 
as described above. If there is no mass concentration directly outside the data field, the second 
method converges faster, while the first approach should be used if the data field does not encom- 
pass most of the mass concentration. Of course, both methods finally approach solutions with 
G x w 1. If magnification information is available, G r should finally approximate unity, which 
provides a useful consistency check of the result. 

3.2 Simulation parameters 

We carry out simulations in which background galaxies are lensed by the mass model shown in 
Fig. |l|. It consists of two softened isothermal spheres, with parameters chosen such that the lens is 
sub-critical. We successfully performed simulations with critical clusters as well, but concentrate 
on the non-critical case to simplify comparisons with the noise-filter reconstructions of Seitz & 
Schneider (1996). The modifications to the practical implementation necessary for critical clusters 
are given in the Appendix. The data field U is a square with side length 5'. We choose a 
number density of 50 galaxies per square arc minute, approximately corresponding to the number 
density which can be achieved with several hours exposure time at a four-meter class telescope in 
good seeing. All reconstructions were performed on a grid with N x — N y = 41. Note that the 
parameters are chosen such that the grid cells have about the same size as the mean separation 
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Figure 1: A mass model consisting of two softened isothermal spheres with core sizes 0'.6 and 0'.8, 
and central surface mass densities k c = 0.4 and 0.5, respectively. Contours are spaced by 0.05, 
and the heavy contour follows k = 0.5 



between galaxies. The intrinsic ellipticity distribution was chosen to be approximately Gaussian, 

= ^[l-cxpM/i^)] ^ I**' 2 ^ ' ™ 

with R = 0.3, yielding M 2 « 0.09, fn w 0.87, and ^ 2 « 0.55 in (jy]) and (Q). The distribution 
p s of relative source sizes is characterized by ay = 0.5 (see Bartelmann & Narayan 1995 for a 
discussion of this choice) . The "observed" ellipticities and relative source sizes are calculated from 
realizations of the source distributions, using the gravitational lens equations. 

In order to assess the expected deviation of G x and G r from unity, we plot in Fig. || the 
quantities 



G x (</3) 

G r (< P) 




§k\<p 



/3-A/3<|fl fc |</3 
l2 \ 



[lnr(fc) - (lnr) (k)Y 

a * I \S k \<p 



GAP) := /^k)-^rHk)?\ . (21) 

\ r / l3-A/3<\e k \<l3 

They are the contributions to G x and G r from galaxies closer than (3 to the center of the data 
field, or within rings of width A/3 = 0'.25 around the center of the data field. Fig. || shows that 
the quantities (^l|) vary considerably between realizations, owing to the broad distribution of 
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Figure 2: Panels (a) through (c) show 
the chi-square contributions of the ob- 
servables, viz. ellipticity and image 
size, within circles of radius (3, and 
within rings of outer radius [3 and 
width A/3 = 0'.25, for three different 
realizations of the source population. 
For this figure, the true mass distri- 
bution of the lens was taken, so that 
these curves merely show the noise 
properties of the source population. 
Heavy solid curves: G x (< (3); light 
solid curves: G x (f3); heavy dashed 
curves: G r (< (3); light dashed curves: 
G r {(3). Panels (d) and (e) show 
the average of these four quantities, 
obtained from 50 realizations of the 
source population, together with their 
1-CT- variations, shown by the same line 
types 

intrinsic source properties. When averaged over 50 realizations (Figs. ||d,e), their mean values are 
very close to unity, and their l-a variations are a good indicator for the expected values in true 
reconstructions. 

Even when the exact deflection potential is used, the resulting mass distribution k will deviate 
from the true surface mass density because k is calculated from ip with finite differencing. For 
the mass distribution shown in Fig. [l], k from finite differencing deviates from the true n by 
< 0.03 everywhere, with the largest deviations occurring at the two mass peaks. The grid is too 
coarse for a more accurate calculation of the second-order derivatives. Since the deviations are 
sufficiently small (i.e. much smaller than the expected accuracy that we can hope to achieve from 
our reconstruction), we have chosen not to further refine the grid, ft should be noted that the 
method by Bridle et al. (f998) suffers from the same, or worse, inaccuracies, because there the 
shear is calculated from the surface mass density by integrating over a coarse grid. 

Finally, Bridle et al. (f 998) calculated the covariance matrix for the resulting mass distribution, 
which wc do not repeat here. One must take into account, though, that the error estimates of the 
resulting mass reconstruction are strongly correlated, because the shear depends non-locally on k. 

4 Results 

4.1 ML reconstructions without magnification information 

We first neglect magnification information, i.e. we minimize (|l3|). For the mass model in Fig. |l|, 
mass reconstructions for 50 realizations of the galaxy population were performed. For each re- 
alization, f 000 iterations were taken, in each case using the noise-filter reconstruction of Seitz & 
Schneider (1996) as initial potential. The number of iterations was chosen to produce a stable 
result for all realizations. The actually required number can be substantially smaller in individual 
cases. Typically, 1000 iterations take approximately 30 minutes on an IBM 590 workstation. After 
every 20 iterations, the prior was changed to the current mass distribution, smoothed by a Gaus- 
sian of width # sm = 0'.2 for the first 600 iteration steps, and 9 sm = 0'.15 afterwards. A somewhat 
larger smoothing length at the beginning leads to faster convergence, but produces artifacts due 
the finite region over which the integral in @ is performed. Note that 9 sm is of the same size 
as the grid cells. In order to avoid excessive fine-tuning for these simulations, the regularization 
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Figure 3: Panels (a) through (c) show 
the chi-square quantities (|2l| ), to- 
gether with the mass ratios (|2 2|) for the 
reconstructions obtained from three 
different realizations of the galaxy 
population. Heavy solid curves: G x (< 
f3): light solid curves: G x (/3); heavy 
dashed curves: G r (< (3); light dashed 
curves: G r (/3); heavy dash-dotted 
curves: m(< (3); light dash-dotted 
curves: m(/3). In panels (d) through 
(f), the average of these quantities 
from 50 realizations are displayed to- 
gether with their l-a range (shown by 
the same line types). For these sim- 
ulations, the regularization parameter 
was kept constant, i.e. not adjusted for 
each simulation to yield G x « 1. No 
magnification information was used 



parameter r/ was fixed to ?y = 30 for all 50 realizations. Of course r\ could be changed to finally 
achieve G x — 1 for each individual realization. The current choice of 77 was made to achieve 
G x w 1 for all realizations (see Fig. || below). The final 20 iteration steps were performed with 
a value of r\ large enough to make the mass distribution follow the prior very closely, and which 
yields a smoothing of the mass reconstruction on a scale of 9 sm . 

Since no magnification information was used in these simulations, the resulting deflection 
potential (and thus the mass distribution) is determined only up to the mass-sheet transformation 
(|J) . In order to compare the reconstructions with the input model, each mass map was transformed 
such that the total mass inside U agreed with the true total mass. In Fig. ]§[ we show the quantities 
( pl| ) for three different realizations of the galaxy population, together with the ratios 

m(B) ■= MW 
[P) ■ M true (/3) ' 

M truo (< p) 

of the reconstructed mass inside rings and circles, relative to the true mass distribution. As can 
be seen, these mass ratios are always very close to unity, which means that the mass maps are 
reconstructed with high accuracy (up to the mass-sheet degeneracy). The dispersion of m{(3) 
about unity is less than 5%, and the mean of m(/3) over 50 realizations is astonishingly flat. 
There is no indication that the mass at the center or in the outer parts is systematically over- or 
underestimated. The choice of the regularization parameter results in G x being slightly smaller 
than unity on average, though with substantial variation from case to case. Evidently, G x {0) is 
significantly smaller in the inner part than in the outer part, an effect also seen in Fig. ^ where the 
true mass distribution was considered. This is due to the fact that the ellipticity distribution after 
lensing is not really a Gaussian, and the deviation from this assumed functional form becomes 
larger for larger values of the reduced shear, i.e. closer to the center of U. Whereas there is no 
fundamental difficulty in replacing the Gaussian with a more accurate probability distribution, 
the simple form for G x is computationally convenient and seems to be sufficiently accurate for the 
mass reconstructions, as seen from the dash-dotted curves in Fig. § The deviation of G r (/3) from 
unity can be substantial in individual reconstructions, but the mean over all realizations is very 
close to unity. 

At the edge of the data field, two systematic effects become visible in the mean quantities 
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plotted in Figs. [|]d-f (and also by looking at the 2-d distribution over IX): The value of G X (JS) 
shows a small but significant decrease, and k is slightly too large near the boundary. The first of 
these effects can be understood by considering the number of galaxies for which the shear estimate 
is affected if the value of ip is changed at a grid point. If that grid point is located in the inner 
part of U, the estimates of 7 for galaxies within the neighboring 16 grid cells are affected. This 
number decreases for points near the boundary of U, so that ip there is less constrained by the 
measured image ellipticities. This implies that at the boundary, it is easier to "fit the noise" 
caused by the intrinsic ellipticity distribution. The slightly too large k near the boundary is due 
to the prior. The prior is obtained from local averaging of the current mass distribution. If the 
mass distribution decreases outwards, the local mean value which can only be taken from the 
grid points within IA will be slightly too large at the boundary, which explains why the method 
presented here is slightly biasing the mass map at the boundary. In actual applications, a strip of 
width 9 snl can be ignored in the analysis of the mass distribution if this bias is a worry. However, 
its amplitude is very small, and it can probably be safely ignored in most situations compared 
to the stochastic errors. Alternatively, one can use a mild extrapolation of the smoothed mass 
distribution to obtain an estimate of the smoothed k values on the boundary less affected by this 
bias. In the case of our mass model, a simple fix (instead of a more elaborated extrapolation) can 
be obtained by decreasing k by An = 0.0015 on all boundary points, which practically eliminates 
the bias. We use this simple fix to the simulations discussed in the next subsection. We further 
point out that the amplitude of this mass bias can also be checked with real data, by generating 
artificial data sets from the reconstructed mass distribution and by performing reconstructions for 
those in the same way in which the original mass reconstruction was obtained. 

0.008 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
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Figure 4: Power spectra k 2 P(k) of the difference between the values of K{9) — ln[l — k(9)] from 
reconstructed mass distributions and the true mass distribution (see Seitz & Schneider 1996 for 
a definition of P{k) and its practical calculation). The solid curve corresponds to if ml, the 
short-dashed curve to Knf- In addition, we plot the power spectra of the difference between the 
ensemble-averaged mass maps, both from the ML method (long-dashed curve) and the noise-filter 
reconstruction (dotted curve). For comparison, the dash-dotted curve shows the power spectrum 
of the difference obtained from the finite-differencing error, mentioned at the end of Sect. 3 

The mass-sheet transformation (^) allows to determine the quantity if (9) = ln[l — k(9)] up 
to an additive constant. As in Seitz & Schneider (1996, 1998; SK), we perform a power-spectrum 
analysis of the difference between if ml and if true, i-e. between the if from the ML reconstructions 
and from the true mass distribution. Note that this power spectrum for k 7^ is independent 
of the transformation (^). In Fig. ^, we compare the power spectrum obtained from our 50 
realizations to the power spectrum obtained from the difference between i^NF and if true, the 
former one corresponding to the noise-filter mass reconstruction of Seitz & Schneider (1996), using 
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the same realizations of the galaxy population. As can easily be seen, the ML reconstruction yields 
substantially lower values for P(k) than the noise-filter reconstruction. 

For illuminating the main reason for this difference, we also calculate the power spectra of the 
difference between (Kml), (-Knf) and -Ktrue, the former two being the mean reconstructed mass 
distribution averaged over all 50 realizations. The difference between these two power spectra is 
very large, with the power of (-Knf) — true being much larger than for the ML reconstructions. 
Almost the entire difference comes from the inability of the noise-filter method to resolve the core 
of the two mass peaks in the model appropriately because of the unavoidable smoothing. Note that 
the smoothing scale AO in the noise-filter reconstruction is much larger than the smoothing scale 
6* sm for the prior in the ML method, the latter being of the same order as the size of a grid cell, or 
the mean separation between two background galaxy images, which in any case is smallest scale 
one can hope to resolve. In fact, it seems that the power spectrum of K^f ~ Ktrae lS approximately 
the sum of the power spectra of Kml — K true and (.Knf) — K tluc . In other words, the error of the 
noise-filter reconstruction is due to a systematic error (from smoothing, although it is adaptive 
smoothing as in Seitz & Schneider 1996), and a stochastic error which is approximately the same 
as that for the ML reconstruction. 

The power spectrum of the error resulting from finite differencing is shown as the dash-dotted 
curve in Fig. |] and is seen to be about the same as that of (if ml) — K tIUC - I.e., the largest part of 
the error of the mean mass map in the ML reconstruction comes from finite differencing; only for 
k < 3 is the latter error larger than that from finite differencing. This little power excess is due 
to the aforementioned slight bias at the boundary of U. 

To conclude this subsection, we have shown that the ML method yields more accurate mass 
reconstructions than the noise-filter method, whereby the largest difference is due to the better 
angular resolution of the ML method in those parts of the data field where the data require more 
local structure, whereas the noise-filter method - like all direct methods developed so far - uses 
a fixed smoothing scale, or at best an ad-hoc choice for the change of smoothing scale with the 
signal, such as in Seitz et al. (1996, and references therein). The ML mass distributions are 
slightly biased at the boundary of the data field, and the amplitude of this bias depends on the 
mass distribution. 

4.2 ML reconstructions with magnification information 

For the 50 realizations of the background galaxies, we also performed ML mass reconstructions 
including magnification information, i.e. minimizing (|l6|). Each reconstruction proceeded as de- 
scribed in the beginning of Sect. 4.1, except that the regularization parameter was chosen to be 
t] = 60 instead of 30. This is necessary to give the regularization the same weight as before, 
because the chi-square part of E is now effectively doubled. The number of iteration steps was 
slightly increased to 1200. Anticipating that the additional magnification information would allow 
a slightly higher angular resolution of the mass map, we somewhat reduced the smoothing scale 
for the prior to 8 sm — O'.l. 

In analogy to Fig. |3j, we plot in Fig. ||the quantities ( ^lf ) and (|3|) for three different realizations 
of the galaxy population, their average over 50 realizations, and their dispersion. Comparing 
the two figures, it seems that including magnification information does not change the curves 
significantly, apart from the fact that the mass-sheet transformation is now obsolete. Magnification 
information therefore mainly affects the normalization of the mass map, but does not provide 
significant information on its shape. 

To further illustrate this point, we plot in Fig. || the power spectra of the difference between 
reconstructed mass maps Kml and the model mass distribution -Ktruc, similar to Fig. || The dif- 
ference in the power spectra between reconstructions with and without magnification information 
is much smaller than their difference to the noise-filter mass reconstructions, at least for small 
values of k < 7. For larger k, the difference becomes somewhat larger, but this may simply be due 
to the different choice of the regularization parameter which affects the small-scale noise. In our 
model, these small scales are only noise because the mass model is relatively smooth. The power 
spectra of the ensemble-averaged mass maps show that the spectrum for simulations with magni- 
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Figure 5: In analogy to Fig. ||, wc 
plot in panels (a) through (c) the 
chi-square contributions (pl|), together 
with the mass ratios (p2[), for the ML- 
reconstructions obtained from three 
different realizations of the galaxy 
population, now using magnification 
information. Heavy solid curves: 
G x (< /?), light solid curves: G x (/3); 
heavy dashed curves: G r {< /3); light 
dashed curves: G r (f3); heavy dash- 
dotted curves: m(< (3); light dash- 
dotted curves: m(/3). In panels (d) 
through (f), the average of these quan- 
tities over 50 realizations is displayed 
together with their \-a range. For 
these simulations, the regularization 
parameter was kept constant, i.e. not 
adjusted for each simulation to achieve 
G v « 1 



fication information is somewhat lower. This is probably due to the fact that the magnification 
information leads to a slightly better angular resolution in the core of the two mass concentrations. 

Our simulation results therefore allow to conclude that magnification information does not 
contribute much information on the shape of the mass distribution compared to shear information. 
For estimating the relative accuracy of the mass determination, one can then assume that the mass 
map is known from the galaxy shapes up to the mass-sheet transformation. Let k be the true 
mass distribution and k(A) the distribution transformed with (^). The magnification information 
can then be used to determine A such as to best fit the galaxy sizes, i.e. to minimize G r . Noting 
that (||) transforms magnifications like /Lt — > ^i/A 2 , the best-fitting A is 



i nX = -±y 

8 k=l 



lnr(fc)- n(0 k ) 



1/2 



(23) 



Since the distribution of the In r is Gaussian, the values of In A for different realizations of the 
galaxy population follow a Gaussian probability distribution with zero mean and dispersion <j\ = 
ay/ \J l ~Ng. Whereas the distribution of In A is symmetric about zero, the mean of A is not unity. 
The ratio 



M(A) (k(A)> 1 / 1 

m(A)= ^r = ^r = ^ + H 1 "w 

of the estimated mass to the true mass has a median of unity, but a slightly smaller mean, 



(24) 



1 - 



dlnAp(ln A) m(ln A) 



- 1 



2N e 



(25) 



The deviation of (m) from unity is very small: For N g = 1250 and ay 
with (k) ~ 0.2, we find (m) « 0.9996. Similarly, the dispersion of m is 



0.5, and a mass model 



1 



- 1 



ay 



(26) 
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Figure 6: Power spectra k 2 P(k) of the difference between the reconstructed mass distributions 
K{9) = ln[l — k(0)] and the true mass distribution. The solid and long-dash-dotted curves show 
Kml obtained without and with magnification information, respectively. For comparison, the 
short-dashed curve shows -FCnf- In addition, we plot the power spectra of the difference between 
the ensemble-averaged mass distributions, both from the ML method without (long-dashed curve) 
and with (dotted curve) magnification information 

where we have used ay / -J N g <C 1 in the approximations. For the same parameters as before, 
Cm = 0.057. We therefore expect that our mass model can be reconstructed with an accuracy of 
about 6% if the shape of the mass distribution is sufficiently constrained by image shapes. Fig. |^d 
shows that this estimate in very good agreement with the mass estimate obtained from the ML 
reconstruction using the magnification information. 

To elaborate on this point, we plot in Fig. the ratios of four different mass estimates to the true 
mass. These are: the mass from ML reconstructions (i) with and (ii) without magnification infor- 
mation, but transformed according to ( p3| ) to best fit the source sizes; (iii) mass-sheet-transformed 
noise-filter reconstructions; and (iv) the mass-sheet-transformed true density distribution. Obvi- 
ously, the scatter of the mass estimates across the four different methods is very much smaller than 
that caused by different realizations of the source population. As much as 50 realizations allow to 
conclude, all four methods yield an unbiased estimate of the total mass and have approximately 
the same dispersion, close to the theoretical expectation discussed above. 

We investigated how the quality of mass estimates is affected by the spatial resolution. The 
noise-filter reconstructions have a lower resolution than the ML reconstructions, although our 
adaptive smoothing makes them satisfactorily resolve the mass peaks. We performed noise-filter 
reconstructions with fixed smoothing scales of 0'.3 and 0'.5. In those cases, the peaks are fairly much 
smoothed out. Applying then ( p3| ) to match image sizes, we found that the total mass within U 
is systematically overestimated, by ~ 1% and ~ 3%, respectively. This is because the smoothing 
of the map renders the model magnification too low to match the observed high magnifications 
close to mass peaks, for which the invariance transformation tends to compensate with a slightly 
higher mass. This bias also affects the ML reconstruction of BNSS because of missing resolution 
in the cluster center. However, it is sufficiently small to be not of much practical importance for 
reconstructions which can be expected from real data in the near future, but it should be kept in 
mind if a large sample of clusters is investigated statistically. 
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Figure 7: For 50 realizations of the galaxy population, we estimate the total mass inside hi and 
compare it to the true mass. Defining to := M/Mtmc, we plot (i) the mass ratio toml obtained 
from ML reconstructions including magnification information, (ii) the ratio rn' Mlj obtained from 
ML reconstructions without magnification information, but mass-sheet transformed according to 
( p3|) , and (iii) the mass ratio TOnf obtained from noise-filter reconstructions, transformed according 
to (^3|), against (iv) the ratio m' tmc obtained from the true shape of the mass distribution, also 
invariance transformed according to (p3|). The means and dispersions of the four mass ratios are: 
(i) (m ML ) = 0.994, cr m = 0.051; (ii) (m' ML ) = 1.000, a m = 0.052; (iii) (m NF ) = 1.002, cr m = 0.051; 
and (iv) (m' tIUC ) = 1.001, a m — 0.054. Within the accuracy allowed by 50 realizations, they all 
agree with (to) = 1 and a m w 0.052, very much in agreement with the estimate ( |26| ) 

5 Discussion and conclusions 

We presented a new method for reconstructing projected mass distributions of galaxy clusters. 
The method uses image distortions of background galaxies and their size as a function of surface 
brightness. Our entropy-regularized ML method (Seitz 1997) is a further development of previously 
published inverse methods for the mass reconstruction. In particular, we describe the lens by its 
deflection potential ip as suggested by BNSS. This is of major importance, for two reasons. First, 
if the surface mass density k on a finite field U is used to describe the lens, the shear on U is 
incompletely specified by the model because the mass distribution outside U can contribute to the 
shear. Second, the shear at the position of any galaxy depends only locally on if>, which allows a 
much faster minimization algorithm for a given number of grid points. 

We regularize the method by an entropy term as suggested by Bridle et al. (1998), but ad- 
ditionally adapt the prior to the current model of the mass distribution. This 'moving prior' 
(Lucy 1994) allows a considerably higher resolution of mass peaks. The spatial resolution of the 
entropy-regularized ML method adapts itself to the strength of the lensing signal, producing mass 
distributions which are as smooth as possible, and as structured as the data require. In that re- 
spect, our method differs from that of BNSS and SK. We showed that the ML method is superior 
to the noise-filter method (Seitz & Schneider 1996) which was the most accurate of the presently 
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known direct inversion methods (Seitz & Schneider 1996, 1998; SK; Lombardi & Bertin 1998). 

Obviously, the method described here is not restricted to rectangular data fields, but can easily 
be adapted to any geometry of U by covering 11 with quadratic grid cells, and adding a boundary 
of grid points for ip - the rest is only a matter of labeling. Furthermore, we note that observational 
errors can be incorporated into the likelihood function. For example, if the measurement error of 
the ellipticity x is fobs, one can replace cr£ in (|l3| and by ct£ + al hs . 

In contrast to the direct inversion methods, all of which are variants and generalizations of 
the original Kaiser & Squires (1993) method, the inverse methods allow to include additional 
information on top of the shear measured through image ellipticities. We demonstrated this 
here by adding magnification information derived from image sizes at given surface brightness, 
as discussed by Bartelmann & Narayan (1995). However, we could equally well use the change 
of number counts due to magnification bias (Broadhurst et al. 1995) as an additional constraint. 
In that case, if the number counts of a certain (e.g. color-selected) galaxy population have a 
cumulative slope of —f3, the expected number density of background galaxies at a position 

is n(9) = n /i(9) , where no are the counts at the same flux limit in the absence of lensing. 
Assuming that galaxies are intrinsically randomly distributed, the probability of having N galaxies 
within U is a Poisson distribution Pn((N)) with 

(AO = n f d 2 9 fx(d) . (27) 
Ju 

Consequently, the likelihood function could be augmented by a factor 

C»=P N {{N))\[U6 k ) . (28) 

If galaxy clustering is important, the likelihood cannot be written as a simple product over 
individual galaxies, but the joint probability distributions must be taken into account. The contri- 
bution of clustering effects to the likelihood function is somewhat uncertain, because an approxi- 
mate expression has to be used due to lack of knowledge on the iV-point correlation functions (see 
Broadhurst et al. 1995 for further discussion). 

Perhaps the most promising generalization of our method is the inclusion of strong lensing 
constraints. Since giant arcs and multiple images of background galaxies provide (nearly) exact 
constraints on the lens mass distribution, it is highly desirable to include them into a mass re- 
construction. The obvious way to do this would be to augment the function E by a term which 
measures the degree to which multiple images of the same source are mapped back to the same 
position in the source plane. In addition, the surface brightness profile of multiple images of ex- 
tended sources can be incorporated, e.g. in a similar manner as the spatially resolved multiple arc 
in the cluster C10024+16 (Colley at al. 1996). 

In some of the observed clusters, the lensing effects of individual galaxies are visible, in partic- 
ular through deformations of giant arcs. Some of the most prominent examples are the triple arc 
in 0024+16 (Kassiola et al. 1992), the multiple arc systems in A 2218 (Kneib et al. 1996), and the 
distortion of the curvature in the arc of the galaxy cB58 in MS1512+36 (Seitz et al. 1998). But 
even weaker lensing effects of individual (cluster) galaxies can be detected using a combination of 
cluster mass reconstruction and galaxy-galaxy lensing techniques. By adding two free parameters 
to the lens model, such as the mass-to-light ratio of cluster galaxies and their characteristic spatial 
extent, the size of halos of cluster galaxies can be investigated (Natarayan et al. 1997; Geiger & 
Schneider 1998). 
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A Notes on the implementation 

We outline here how the function E - see ( |l3| ) and (|l^) - can be calculated from the deflection 
potential {ip} at the grid points. We first note that E can be easily calculated in terms of the 
shear 7 and the surface mass density k at the position of the galaxies 9 k , and in terms of k at the 
grid points, which enters into the regularization term 1Z. 

Let ip a label the deflection potential at the A^im grid points where ip is varied. If no mag- 
nification information is used, 7V dim = (N x + 2)(N V + 2) - 8, and N dim = (N x + 2)(N y + 2) - 7 
otherwise. The components of 7 and K at the N x N y grid points are obtained by second-order 
finite differencing, according to (|I]) and For calculating 72 at the grid corners, we used the 
finite-difference formula (25.3.27) from Abramowitz & Stegun (1984) so that ip at the corner points 
of the extended grid does not enter E. This can be important because the corner points do not 
enter the entropy term, so that they are not regularized. If no galaxy happens to lie in the corner 
grid cells, these points would be unconstrained. Finite differencing being a linear operation, we 
can write 

K„ = D^ a + d M , 1 < n < N x N y . (29) 

Summation over a, 1 < a < N^i m , is implied, and the array accounts for the contribution to 
from the grid points where ip is held constant. For each galaxy position 6k, the values of k 
and 7 are obtained from bilinear interpolation from the closest grid points. Since interpolation is 
again a linear operation, we can write 

7l(0fc) = A ka ip a + a k , 
l2{0 k ) = B ka ip a + b k , 

K(6k) = C ka ip a + c fe , (30) 



for 1 < k < N g . The arrays a, 6, and c again arise from the contributions of the grid points where 
ip is held fixed. If magnification information is used, ip is kept fixed at three grid points, where 
ip — Q can be chosen without loss of generality; then, a k = b k = c k = = d M . If magnification 
information is neglected, at least one of the four constant values of ip has to differ from zero, and a, 
b, c, and d are not identically zero, although most of their elements vanish. The same holds for the 
elements of the matrices A, B, C, and D. Whereas the dimension of these matrices is large, only 
about 12Af g elements are different from zero. Thus, we store these matrices using the row-indexed 
sparse storage mode, as described in Press et al. (1992, Sect. 2.7). Note that these matrices and 
'vectors' have to be calculated only once, since they depend only on the number of grid points 
and the location of the galaxies. The sparseness of the matrices is one of the advantages to work 
in terms of the deflection potential rather than the surface mass density. 
For the calculation of dE/dip a , we note that for any function f k at 9 k , 



df k df k 



dipa 7l (0 fc ) 
For the regularization term, we find 



-A kn 



dfk 

72(4) 



-Bi. 



dfk 

K(9 k ) 



■C k 



(31) 
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dlp a 



N x N y 



N X Ny , , 

fj,= l V ^ 



(32) 



Using these relations, the function E and its derivative with respect to the ip a are easily calculated. 

An additional complication arises for a critical lens. As can be seen from the explicit form of 
the terms of E, the likelihood function attains singularities if a galaxy image happens to lie on 
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a critical curve where \8\ = 1 and yT 1 = 0. This means that the minimization procedure cannot 
modify the model such that critical curves move over galaxy positions. Since the initial guess of 
tjj will certainly not be accurate enough that such crossings can be avoided, we have to modify 
the function E appropriately. The necessary modification of E can be achieved by replacing 

by [m~ 2 + f2 ] > an d &x by °x + V) where e and rj are two small quantities. These replacements 
leave E finite if a galaxy image is placed on a critical curve. The minimization then proceeds by 
setting e and rj to about 0.1 at the beginning of the minimization, and then slowly decreasing them 
in later iteration steps. This leads to convergence without additional problems. In addition, as 
was true for the direct inversions, by considering a broad redshift distribution of the background 
galaxies (Seitz & Schneider 1997), the singularities connected with critical curves are avoided 
(Geiger & Schneider 1998). 
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